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The theoretical effects of phase separation on encounter-limited charge carrier recombination in 
organic semiconductor blends are investigated using kinetic Monte Carlo (KMC) simulations of 
pump-probe experiments. Using model bulk heterojunction morphologies, the dependence of the 
recombination rate on domain size and charge carrier mobility are quantified. Unifying competing 
models and simulation results, we show that the mobility dependence of the recombination rate can 
be described using the power mean of the electron and hole mobilities with a domain size dependent 
exponent. Additionally, for domain sizes typical of organic photovoltaic devices, we find that phase 
separation reduces the recombination rate by less than one order of magnitude compared to the 
Langevin model and that the mobility dependence can be approximated by the geometric mean. 
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Despite vast scientific investigations on organic elec¬ 
tronic devices over the last two decades, significant gaps 
in fundamental understanding exist in key areas. A de¬ 
tailed description of charge carrier recombination, im¬ 
portant for designing organic photovoltaics (OPVs), or¬ 
ganic light-emitting diodes, and organic photodiodes, is 
a work in progress. To continue improving these devices, 
the dominant factors controlling recombination processes 
must be understood further so that they can be carefully 
controlled in energy efficient devices. 

Bimolecular charge recombination in organic semicon¬ 
ductors is most commonly described as a second-order 
reaction following the Langevin model, [1] 

R = fc L np, (1) 

where /cl is the Langevin recombination coefficient and 
n and p are the concentrations of electrons and holes, 
respectively, k l is derived by assuming an encounter- 
limited reaction in which the time it takes for an elec¬ 
tron and hole to come together due to their Coulomb 
attraction is rate limiting. As a result, 

g 

kh — (Pe T ph)] (2) 

where e is the elementary charge, e is the dielectric con¬ 
stant, eo is the vacuum permittivity, and p e and ph are 
the electron and hole mobilities, respectively. 

The Langevin model also assumes a spatially and ener¬ 
getically homogeneous and isotropic system with no in¬ 
ternal electric held, which is not strictly valid in most 
organic semiconducting devices. Organic semiconduc¬ 
tors are well-characterized as having varying degrees of 
energetic and spatial disorder that can have a major im¬ 
pact on charge transport properties as highlighted by the 
commonly used Gaussian disorder model (GDM).[2] In 
addition, devices may operate with a significant internal 


electric held. Investigating these issues, a number of ki¬ 
netic Monte Carlo (KMC) simulation studies have identi¬ 
fied conditions where deviations from the Langevin model 
occur.[3-7] However, van der Holst et al. concluded that 
the Langevin model still works well in an isotropic system 
with the GDM at low electric helds as long as accurate 
mobility values are used. [8] In agreement, measurements 
on a number of neat organic semiconducting materials 
have been consistent with the Langevin model. [9-12] 
While the Langevin model may work well for neat ma¬ 
terials, many devices utilize phase separated blends. [13] 
In such blends, electrons and holes are relegated to sepa¬ 
rate phases and are only able to undergo recombination 
at the phase boundaries. The resulting spatial limitations 
on charge carrier motion and recombination locations 
are expected to alter the recombination kinetics. Bulk 
hetero junction (BHJ) OPVs with domain sizes ranging 
from approximately 10 to 50 nrn have often exhibited 
two major deviations from the Langevin model. First, 
super-second order recombination kinetics has been mea¬ 
sured in several blend systems.[14-20] Often attributed 
to charge traps, the reasons for this behavior are still 
under debate.[12, 19, 21-26] 

The second major deviation, commonly observed in 
P3HT:PCBM films, is a recombination rate that is sev¬ 
eral orders of magnitude less than predicted by the 
Langevin model.[16, 27-29] As a result, Pivrikas et al. 
proposed that a reduced recombination rate is an inher¬ 
ent property of BHJ blends,[27] and Koster et al. created 
the minimum mobility model, arguing that the recombi¬ 
nation rate in a phase separated system should be limited 
by the mobility of the slowest carrier, [30] 

g 

k m in = —min(/i e , /i^). (3) 
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However, even though a significant reduction was also 
found in several other blends,[31-36] some blend systems 
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exhibit recombination rates much closer to the Langevin 
model. [31, 32, 37-39] It is becoming increasingly clear 
that a greatly reduced recombination rate is not an in¬ 
herent property of phase separated blends but is a prop¬ 
erty that is dependent on a number of factors that are 
still under debate.[12, 23, 33, 40-44] 

Despite clear deviations in many cases, the Langevin 
model is often used to explain bimolecular recombina¬ 
tion in phase separated blends because a more complete 
model is still missing. An important first step towards a 
more complete model is understanding the fundamental 
effect of phase separation on encounter-limited bimolec¬ 
ular recombination. Previously, Groves and Greenham 
used KMC simulations to show that the recombination 
rate in a simple phase separated system lies somewhere 
between the Langevin model and the minimum mobility 
model with a weak dependence on domain size. [7] In ad¬ 
dition, some experiments have also indicated a relatively 
weak domain size dependence. [34] However, other KMC 
simulations[45] and experiments[36] have indicated that 
the domain size could have a larger impact. 

In this letter, we determine the dependence of the re¬ 
combination rate on the domain size and the electron 
and hole mobilities. We find that with a very small do¬ 
main size, the Langevin model still holds, but for larger 
domains, clear deviations are present. Unifying the com¬ 
peting Langevin and minimum mobility models with our 
simulation results, we show that the mobility dependence 
can be described using the power mean of the mobilities 
with a domain size dependent exponent. Additionally, 
for domain sizes typical of OPVs, we find that phase 
separation reduces the recombination rate by less than 
one order of magnitude compared to the Langevin model 
and that the mobility dependence can be approximated 
by the geometric mean. 

Due to the complex geometry of phase separated sys¬ 
tems, analytical derivation of the recombination rate as 
a function of the domain size is extremely difficult. As 
a result, KMC simulations were performed to reach a 
numerical solution. The simulations were configured to 
simulate pump-probe experiments on BHJ films without 
electrodes. Using the Ising_OPV software tool, a simple 
Ising phase separation model was implemented to create 
morphologies with varying domain size. [46, 47] To verify 
the presence of each recombination regime, the domain 
size was varied from 5 to 55 nm. For each domain size 
(d) tested, 100 morphologies were independently gener¬ 
ated to form a morphology set. The resulting morpholo¬ 
gies consist of two pure phases with equivalent average 
domain sizes in a bicontinuous network configuration. 

The model morphologies were then implemented into a 
three-dimensional lattice with a lattice constant of 1 nm, 
and sites were assigned energies from an uncorrelated 
Gaussian DOS with an energetic disorder (a) of 75 meV. 
Two-dimensional periodic boundary conditions were used 
to simulate a thin film. To start the simulation, excitons 


were created with uniform probability throughout the 
lattice with a Gaussian excitation pulse having a pulse 
width of 100 ps and an intensity corresponding to an ini¬ 
tial exciton concentration of 5 x 10 17 cm -3 . However, 
for 55 nm domains, the initial exciton concentration was 
set to 3.1 x 10 17 cm -3 due to computational limitations. 
This change had no impact on the recombination behav¬ 
ior of interest in this study. The complexities of charge 
separation were bypassed to create free charge carriers di¬ 
rectly from excitons. To do this, electron-hole pairs were 
created across the interface with a separation distance of 
30 nm by restricting exciton creation to within 30 nm of 
an interface and executing an ultrafast long-range charge 
transfer event. The charges in the lattice then under¬ 
went hopping transport using the Miller-Abrahams (MA) 
model, [48] and coulomb interactions were included be¬ 
tween charges within a cutoff radius of 35 nm. Electron 
hopping was restricted to acceptor domains and hole hop¬ 
ping was restricted to donor domains. Charge recombi¬ 
nation was also implemented using the MA model with a 
recombination prefactor (i?o,rec) that was held constant 
at a large value of 10 15 s -1 to ensure that recombination 
dominated over redissociation. 

For each simulation, 24 morphologies were randomly 
selected from the appropriate morphology set and 4 
random configurations of energetic disorder were imple¬ 
mented for each, and the results of the 96 runs were av¬ 
eraged. During each simulation, the hole concentration 
was logged as a function of simulated time. With the 
lattice sizes used, the carrier concentration could be re¬ 
solved over two and half orders of magnitude, covering 
a range typical for steady state illumination intensities 
from 0.1 to 10 suns. Assuming second-order recombina¬ 
tion kinetics and n = p, the numerical derivative of the 
hole concentration as a function of time ( t ) was then used 
to calculate the simulated time-dependent recombination 
coefficient: 


dp(t) 

fcsim(t) = ~p{t P ^ 

In addition, the displacement of each carrier from its 
initial position was recorded over its lifetime, and the 
numerical derivative of the average mean squared dis¬ 
placement over time for all carriers was used to calculate 
the average time-dependent diffusion coefficient. Due 
to the thin film geometry, the two-dimensional diffusion 
equation was found to be most appropriate. Using the 
Einstein relation, which remains valid at zero-field[4, 49] 
when recombination removes deeply trapped carriers, [50] 
the average time-dependent zero-field mobility of each 
carrier type was determined, 


V(t) 


e d{r(t) 2 ) 
4/cbT dt 


( 5 ) 


More information about the morphology generation, 
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KMC simulation parameters, and data analysis is shown 
in the Supplementary Information. [51] 

To determine the mobility dependence, the relative 
magnitudes of the electron and hole mobilities were tuned 
by varying the hole hopping prefactor (I?o,h) from 10 11 
to 10 15 s _1 while holding the electron hopping prefactor 
(Ro, e ) constant at 10 13 s -1 . Figure 1 shows how the hole 
concentration ( p ) decays over time for differing magni¬ 
tudes of Ro,h- By the time the hole concentration reached 
10 16 cm -3 (a typical concentration for 1 sun steady state 
illumination), all transients showed steady second-order 
decay. The time point where p = 10 16 cm -3 was used as 
the comparison point between different simulations. At 
this time point, the simulated recombination coefficient 
(fcsim), the mobilities (p, e and the Langevin recombi¬ 
nation coefficient ( k l), and the minimum mobility model 
recombination coefficient (fc m ; n ) were calculated. 



Time (s) 

FIG. 1. (a) Hole concentration transients and (b) time- 

dependent mobilities for d = 15 nm with varying hole hopping 
rates. 

Figure 2 shows how the mobility dependence of the 
recombination coefficient changes with different magni¬ 
tudes of domain size. With a very small domain size of 
5 nm, the recombination rate approaches the Langevin 
model, and as the domain size increases, the recombi¬ 
nation rate deviates from the Langevin model. In the 
intermediate regime, in agreement with the results of 
Groves and Greenham,[7] the simulated recombination 
rate was not found to be proportional to the sum of mo¬ 


bilities nor the minimum mobility. Instead, a new trend 
not captured by any of the current models is observed, 
in which fc s j m is approximately proportional to the geo¬ 
metric mean of the electron and hole mobilities. For the 
largest domain sizes tested, the mobility dependence con¬ 
tinues to change as observed by the downward curvature 
beginning to form for d=45 and 55 nm. Here, the mo¬ 
bility dependence appears to be approaching the shape 
predicted by the minimum mobility model. However, the 
bias toward the minimum mobility appears much weaker 
than predicted in the minimum mobility model. Another 
type of mean which favors the smaller value is the har¬ 
monic mean, and the simulated recombination behavior 
appears to be closer to the harmonic mean. The same 
general behavior was also observed when using different 
magnitudes of energetic disorder (cr) and when perform¬ 
ing the analysis at carrier concentrations from 5 x 10 15 
to 2 x 10 16 cm -3 . 
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FIG. 2. The effect of the domain size on the mobility depen¬ 
dence of the simulated recombination rate coefficient com¬ 
pared to the Langevin, harmonic mean, and minimum mobil¬ 
ity models. 

To understand this behavior, consider a case where 
the electron has a very fast mobility compared to the 
hole. In a neat material, even if the hole is slow, it is 
still accessible to the electron everywhere. In this simple 
case, the resulting recombination rate is proportional to 
the sum of the individual mobilities as described in the 
Langevin model. Now, let us consider a phase separated 
morphology with very small domains, as shown in Fig. 
3a. In this case, with highly interconnected domains and 
high interfacial area, the slow hole is still accessible to 
the electron almost everywhere. When the electron en¬ 
ters the Coulomb capture radius of the hole (depicted 
by the dashed circle), the attractive force becomes very 
strong and recombination occurs quickly. As a result, 
even though the electron and hole are restricted to sepa- 
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rate phases, the resulting recombination rate should still 
follow the Langevin model. 



FIG. 3. Depiction of the two extreme recombination regimes: 
(a) very small domain size and (b) very large domain size. 
One charge is depicted with a dotted circle around it that rep¬ 
resents the Coulomb capture radius ( Rc )• The other charge 
is shown with a green path that illustrates one of the most 
likely pathways to enter the capture sphere. 


However, with very large domains, as shown in Fig. 
3b, most of the slow holes are inaccessible to the elec¬ 
trons, and in order for recombination to occur, the hole 
must first migrate to the interface. As a result, the re¬ 
combination rate should be limited by the mobility of 
the slower hole and should follow the minimum mobility 
model. Given these two well-defined extreme cases, it 
follows logically that for intermediate domain sizes, the 
recombination behavior should be somewhere inbetween. 
In effect, the recombination rate of the slower carriers de¬ 
pends mainly on their distance from the interface. Those 
located right next to the interface will require almost no 
further motion, and their recombination rate will be dom¬ 
inated by the speed of the faster carrier. However, those 
located farther away will start to be limited more by their 
own slow motion, passing through a regime where the 
recombination rate will depend on the speed of both car¬ 
riers and then eventually dominated by the slow carrier 
when very far away. Given this behavior, when changing 
the domain size, the critical change is in the average dis¬ 
tance between the carriers and the interface, which then 
alters the mobility dependence. 

Capturing this behavior and unifying the competing 
analytical models with the simulations results of Groves 
and Greenham[7] and those presented here, all curves in 
Fig. 2 can be fit by the following simple equation (see 
Supplementary Information) ,[51] 

g 

^sim = fl ((fy'ZMg/d) ^h)i (6) 

where fi(d) is a domain size dependent prefactor and 
M g [d){^e, Hk) is the power mean (generalized mean), 




l4+f*h 


1/9 


( 7 ) 


with a domain size dependent exponent, g(d). For a very 
small domain size, lim^o fi(d) = 1 and lim^o g{d) = 1, 
and the Langevin expression is obtained. In addition, 
lim 0 _>o M g is the geometric mean, M_i is the harmonic 
mean, and M_ 00 is the minimum value. 

Fig. 4 shows the fitted values for /i and g as a func¬ 
tion of the domain size. For d = 5 nm, /i and g both 
approach one, meaning that the behavior is very close 
to the Langevin model. For intermediate domain sizes 
of about 10-35 nm, g « 0, meaning that the behavior in 
this regime can be approximated by the geometric mean. 
Then for d > 35 nm, g continues to slowly decrease, 
thereby increasing the impact of the minimum mobility. 
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FIG. 4. Fitted values for the prefactor (/i) and the power 
mean exponent ( g ) as a function of the domain size. 

In addition to the change in the mobility dependence, 
the magnitude of the recombination rate also decreases 
with increasing domain size. It is clear that simple phase 
separation alone causes only a mild reduction of less than 
one order of magnitude when the electron and hole mobil¬ 
ities are within one order of magnitude of each other, as 
is the case in most optimized devices. As a result, blends 
with much larger reductions must have other more dom¬ 
inant contributing factors. These blends are likely to be 
greatly affected by polaron pair redissociation and are 
therefore not in the encounter-limited regime and/or are 
affected by interfacial states different from those in the 
bulk. [12] 

In conclusion, we have quantified the effect of phase 
separation on the bimolecular charge recombination rate 
when in the encounter-limited regime. Most significantly, 
we have shown that the mobility dependence can be de¬ 
scribed by the power mean with a domain size dependent 
exponent. For domain sizes typical of OPVs, the geomet¬ 
ric mean is a very good approximation of the mobility 
dependence. In addition, we clearly demonstrate that 
greatly reduced recombination rates are not an inher¬ 
ent property of phase separated systems. These results 
represent a unification of previous recombination models 
and a major step forward in developing a more complete 
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model for charge carrier recombination in organic semi¬ 
conductor blends. With this knowledge, updated design 
rules for new materials and device architectures can be 
defined. 
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MORPHOLOGY DETAILS 


hopping, 


Morphologies were created using the Ising_OPV soft¬ 
ware tool [1, 2] with a 50:50 blend ratio and an interac¬ 
tion energy of 0.6 kT. Utilizing the smoothing and rescal¬ 
ing methods, [2] morphologies were generated on lattices 
with a final size of 100 by 100 by 100 or larger to pre¬ 
vent lattice confinement effects. [2] For each set of input 
parameters, 100 morphologies were independently gen¬ 
erated to form a morphology set. Eight morphology 
sets (MS1,MS2,MS3,MS4,MS5,MS6,MS7,MS8) were gen¬ 
erated with different domain sizes of approximately 5, 
10, 15, 20, 25, 35, 45 and 55 nm. More detailed infor¬ 
mation on the input parameters and characterization of 
each morphology set is provided in Table I. 


KMC SIMULATION DETAILS 


The KMC simulation methodology used in this study 
was explained in more detail in a previous paper,[3] but 
a short summary of the most important aspects are pro¬ 
vided here. The model morphologies were used to define 
the donor and acceptor sites on a three-dimensional lat¬ 
tice with a lattice constant (a) of 1 nm. Both phases were 
assigned an uncorrelated Gaussian DOS characterized by 
the energetic disorder parameter (er). Periodic boundary 
conditions were used in two directions to simulate a thin 
film. To start the simulation, excitons were created with 
uniform probability throughout the lattice with a Gaus¬ 
sian excitation pulse having a pulse width of 100 ps and 
an intensity corresponding to an initial exciton concen¬ 
tration of 5 x 10 17 cm -3 . 

Exciton diffusion was implemented using the Forster 
resonance energy transfer model, 


Rij,ex h — h > Q (:x t| /b (^E/ij, ex h), (1) 

where i?o, e xh is the exciton hopping prefactor, d.^ is the 
distance between sites, 


/b(AI^) 



A Eij > 0 
A E i:i < 0 ’ 


( 2 ) 


and AEy exh is the change in potential energy for exciton 


A E-ijexh E,. s j n g| ( „ £7 singlet- (3) 

Exciton hopping events were calculated to sites up to 
4 nm away from the starting site. In additon, the exciton 
relaxation time defines the lifetime of the excited state 
and is used to calculate the exciton relaxation rate, 

f^exr = 1 /"Ax - (4) 

where r ex is the exciton lifetime. 

The complexities of charge separation were bypassed to 
create free charge carriers directly from excitons. To do 
this, electron-hole pairs were created across the interface 
with a separation distance of 30 nm by restricting exciton 
creation to within 30 11m of an interface and executing 
an ultrafast long-range charge transfer event. Long range 
charge transfer (exciton dissociation) was implemented 
using the simplified Miller-Abrahams model where charge 
transfer is always energetically favorable, 

Rij,ex d = T^o.exd exp ( ) (5) 

where f?o, ex d is the exciton dissociation prefactor and 7 ex 
is the inverse exciton localization parameter. Exciton 
dissociation events were only calculated for sites between 
30 and 31 nm away from the starting site. 

Charge motion was simulated using the Miller- 
Abrahams model. For electrons, 

Edj.eih = A(|,. exp( 27(d) tljj ),/Vj(AU jen,) (6) 

where f?o, e is the electron hopping prefactor, y c h is the 
charge localization parameter, and 

AEy ie ih = E7LUM0 — Ej ,lumo + A Ec,ij — Fdij , (7) 

where £Vlumo and Eij lumo are the initial and final site 
energies drawn from the density of states distribution, 
A Ec,ij is the change in Coulomb potential that would 
occur for hopping from site i to site j , and F is the elec¬ 
tric field. Analogous expresions are used to calculate the 
hole hopping rate. Coulomb interactions were included 
between charges within a cutoff radius. The cutoff radius 
should be set to a large enough value such that it does 
not affect the results of the simulation. A cutoff radius 
of 35 nm was found to be large enough to not impact the 
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TABLE I. Morphology Set Information 



MSI 

MS2 

MS3 

MS4 

MS5 

MS6 

MS7 

MS8 

Initial lattice dimensions 

100 

50 

34 

25 

23 

23 

41 

28 

Monte Carlo steps 

369 

374 

374 

369 

374 

369 

6000 

1110 

Rescale factor 

N/A 

2 

3 

4 

5 

7 

5 

3,3 

Smoothing threshold 

0.52 

0.52 

0.52 

0.52 

0.52 

0.52 

0.52 

0.52 

Final lattice dimensions 

100 

100 

102 

100 

115 

161 

205 

252 

Domain size, d 

4.92±0.02 

10.1T0.1 

15.1T0.3 

20T0.7 

25±1 

35±1 

45±2 

55±3 

Interfacial area/volume 

0.356T0.001 

0.175±0.001 

0.116T0.002 

0.087T0.002 

0.069T0.002 

0.050T0.001 

0.037T0.001 

0.031T0.001 

Tortuosity 

1.1T0.02 

1.1T0.03 

1.1T0.04 

1.1T0.05 

1.1T0.05 

1.1T0.05 

1.1T0.06 

1.1T0.06 


recombination rate at a charge carrier concentration of 
10 16 cm -3 . The change in Coulomb potential is calcu¬ 
lated, 


AE c ,ij = Ecj — E c ,i, 

where 

N 

Ec,i = V] . q,qk , d ik < 35 nm, 

*=f Ui 4nee ° dik 


( 8 ) 

( 9 ) 


given N nearby electrons and holes. An analogous ex¬ 
pression is used to calculate the Coulomb potential for 
the final state by assuming that the charge of interest 
is positioned on site j. Electron hopping was restricted 
to acceptor domains and hole hopping was restricted to 
donor domains. Charge hopping events were calculated 
for sites up to 3 nm away from the starting site. 

When an electron and a hole come close together, the 
charge recombination event is enabled. Charge recombi¬ 
nation was also implemented using the Miller-Abrahams 
model similar to charge hopping, 


Eij, rec — A.Q.rcc CXp ( ), ( 10 ) 

where i?o,rec is the recombination prefactor, which was 
held constant at a large value of 10 15 s _1 to ensure that 
recombination dominates over redissociation. Charge re¬ 
combination events were also calculated for sites up to 
3 nm away from the starting site. In addition, the selec¬ 
tive recalculation method [3] was used with a recalcula¬ 
tion cutoff radius of 5 nm. A full list of parameters is 
provided in Table II. 


SUPPLEMENTARY RESULTS 


All simulated recombination coefficient data was fit 
with the equation, 

ksim = — fi{d)2M g{d) (g e , g h ). (11) 

where /i(d) is a domain size dependent prefactor and 
M g ^(n e , fXh) is the power mean (generalized mean), 


Afg (@e; M/l) 


m! + /4 


i /g 


( 12 ) 


with a domain size dependent exponent, g(d). Fig. 1 
shows the resulting fits of each data set. The resulting 
fitted parameters fi and g and their uncertainties are 
shown in the main article. These results show that the 
power mean in Eqn. 11 is a very accurate description of 
the mobility dependence obtained in the recombination 
simulations. 
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TABLE II. KMC Simulation Parameters 


Lattice constant, a 

1 nm 

Temperature, T 

300 K 

Dielectric constant, e 

3.5 

Energetic disorder, a 

0.075 eV 

Exciton lifetime, r ex 

500 ps 

Exciton hopping prefactor, Ro.exh 

10 12 s _1 

Exciton localization, y e x 

0.1 nm -1 

Exciton dissociation prefactor, Ro.exd 

10 16 s- 1 

Electron hopping prefactor, i?o,e 

10 13 s' 1 

Hole hopping prefactor, Ro,h 

variable 

Charge localization, y C h 

2 nm -1 

Charge recombination prefactor, Ro.rec 

10 15 s- 1 

Electric field, F 

0 VnU 1 

Coulomb cutoff radius, -Rcutoff 

35 nm 
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Recombination Coefficient, k sim (cm s ) 


3 



FIG. 1. Fitted recombination coefficient data. 



